El presente reporte técnico recopila el desarrollo del trabajo final de la asignatura Analítica Predictiva dictada por la Facultad de Minas en la Universidad Nacional de Colombia sede Medellín.
## Objetivo General.
Predecir los accidentes en la ciudad de Medellín, tomando como insumo los datos abiertos de movilidad, publicados por la Alcaldía de Medellín en el portal GeoMedellin.
## Entrenamiento de un modelo preditivo.
Se requiere construir un modelo de prediccion que pueda pronosticar los accidentes en la ciudad a nivel de barrio por cada clase de accidente a una periodicidad diaria, semanal y mensual
## Agrupamiento de los barrios de Medellín según su accidentalidad.
Además del modelo se exige desarrollar un agrupamiento de los barrios de Medellín de acuerdo con su accidentalidad.
## Presentación.
Los hallazgos y el modelamiento están disponibles para su uso a través de una aplicación web, la cual puede visitarse “aquí”, que cuenta además con video promocional.
Medellín con 2’569.007 (proyección del DANE 2020) habitantes es la segunda ciudad más poblada de Colombia, la cual es caracterizada por ser la ciudad con mejor calidad de vida, según Ciudades Cómo Vamos; y ser la segunda ciudad en aportar más al PIB del país, al rededor de 14,38% a 2018. La fuerte industria de la ciudad se ha expandido a los municipios metropolitanos, constuyendo de esta manera, una zona del País con alto desarrollo económico.
Tomando como punto de partida la importancia de la ciudad, se espera que al ser un municipio agitado por sus condiciones económicas y sociales, por consiguiente, se tengan cifras de accidentalidad altas, por tanto resulte conveniente desarrollar un modelo de pronóstico que permita predecir la cantidad de accidentes y así llevar a cabo planes de acción que busquen mitigar al máximo la accidentalidad en la ciudad.
Como se menciona en la introducción y en el objetivo general, se trabaja con datos abiertos proporcionados por la Alcaldía de Medellín sobre el registro de accidentes en la ciudad y corregimientos aledaños.
Librerias usadas:
## Captura de Datos.
Para dar cumplimiento a los parámetros establecidos, se descargan las series de datos pertenecientes a los años 2014, 2015, 2016, 2017 y 2018, por consiguiente se procede a combinar todos los archivos para obtener un sólo set de datos.
urlfile <- "https://raw.githubusercontent.com/jdgallegoq/analitica_predictiva/master/accidentes3.csv"
accidentes <- read.csv(urlfile, fileEncoding = "ISO-8859-1")
#accidentes <- read.csv("D:/jhoparra/2020-1/Predictiva/analitica_predictiva/accidentes3.csv",fileEncoding = "ISO-8859-1")
El set de datos original contiene las siguientes variables:
colnames(accidentes)
## [1] "OBJECTID" "RADICADO" "FECHA" "HORA"
## [5] "DIA" "PERIODO" "CLASE" "DIRECCION"
## [9] "DIRECCION_ENC" "CBML" "TIPO_GEOCOD" "GRAVEDAD"
## [13] "BARRIO" "COMUNA" "DISENO" "DIA_NOMBRE"
## [17] "MES" "MES_NOMBRE" "X_MAGNAMED" "Y_MAGNAMED"
## [21] "LONGITUD" "LATITUD"
De las cuales las más relevantes para el desarrollo de todo el trabajo son: 1. Fecha. 2. Día. 3. Periodo (Año). 4. Clase 5. Gravedad. 6. Barrio. 7. Comuna. 8. Dia Nombre. 9. Mes. 10. Mes Nombre. 11. Longitud. 12. Latitud.
las clases de accidentes descritas son las siguientes:
unique(accidentes$CLASE)
## [1] "Atropello" "Choque" "Otro"
## [4] "Caida ocupante" "Volcamiento" "Incendio"
## [7] "Choque y Atropello" ""
Por otra parte las gravedades de los accidentes son:
unique(accidentes$GRAVEDAD)
## [1] "MUERTO" "SOLO DAÑOS" "HERIDO"
Tomando como base las clases y las gravedades de accidentes, por año se tienen las siguientes cifras:
table(accidentes$PERIODO, accidentes$CLASE)
##
## Atropello Caida ocupante Choque Choque y Atropello Incendio Otro
## 2014 0 4779 4157 27157 0 8 4521
## 2015 1 4485 3691 28250 0 1 4218
## 2016 5 4167 3680 28631 0 4 4870
## 2017 0 3638 3433 29196 1 4 4718
## 2018 0 3610 3617 28213 0 7 3739
##
## Volcamiento
## 2014 972
## 2015 1435
## 2016 1484
## 2017 1572
## 2018 1174
table(accidentes$PERIODO, accidentes$GRAVEDAD)
##
## HERIDO MUERTO SOLO DAÑOS
## 2014 23077 256 18261
## 2015 23274 250 18557
## 2016 23994 235 18612
## 2017 22865 226 19471
## 2018 20991 297 19072
El siguiente gráfico de barras apiladas permitirá diferenciar de manera visual las frecuencias de cada gravedad por año.
acc_fig1 <- accidentes[, c("PERIODO", "CLASE", "GRAVEDAD")]
tabla <- table(acc_fig1$PERIODO, acc_fig1$GRAVEDAD)
acc_fig1 <- as.data.frame.matrix(tabla)
colnames(acc_fig1) <- c("herido", "muerto", "solo_danos")
fig1 <- plot_ly(acc_fig1, x = ~unique(accidentes$PERIODO),
y = ~herido, type = 'bar', name = 'Herido', text = acc_fig1$herido,
textposition = 'outside',
marker = list(line = list('rgb(105,105.105)')))
fig1 <- fig1 %>% add_trace(y = ~solo_danos, name = 'Solo Daños', text = acc_fig1$solo_danos,
textposition = 'outside',
marker = list(line = list('rgb(105,105.105)')))
fig1 <- fig1 %>% add_trace(y= ~muerto, name = 'Muertos', text = acc_fig1$muerto,
textposition = 'outside',
marker = list(line = list('rgb(105,105.105)')))
fig1 <- fig1 %>% layout(title = "Accidentes por Clase y Año",
yaxis = list(title = 'Accidentes'),
xaxis = list(title = 'Año'),
barmode = 'stack')
fig1
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
Se puede percibir que el año 2016 ha sido el año con mayor accidentalidad entre 2014 y 2018, sin embargo, los accidentes mortales no son percibibles en las barras apiladas dada la diferencia de frencuencias con las otras gravedades; por esta razón se exponen únicamente los accidentes mortales en el siguiente gráfico:
fig2 <- plot_ly(acc_fig1, x = ~unique(accidentes$PERIODO),
y = ~muerto, type = 'bar', name = 'Accidentes Fatales',
text = acc_fig1$muerto, textposition = 'outside',
marker = list(line = list('rgb(105,105.105)')))
fig2 <- fig2 %>% layout(title = "Accidentes Fatales por año",
yaxis = list(title = 'Acc'),
xaxis = list(title = 'Año'),
barmode = 'stack')
fig2
Aunque el gráfico de barras apiladas mostraba al año 2014 como el año con menos accidentalidad, el año 2014 fue año con mayor cantidad de accidentes fatales.
Después de haber graficado las frecuencias de accidentes por año, la siguiente serie de tiempo expone la evolución de los accidentes por mes, comparando por lineas cada año de resgistros:
acc_fig3 <- as.data.frame(as.Date(accidentes[, "FECHA"], format = "%Y-%m-%d"))
colnames(acc_fig3) <- "FECHA"
acc_fig3$MES <- strftime(acc_fig3$FECHA, format = "%m")
acc_fig3$ANO <- strftime(acc_fig3$FECHA, format = "%Y")
acc_fig3 <- as.data.frame.matrix(table(mes = acc_fig3$MES, ano = acc_fig3$ANO))
acc_fig3$MES <- row.names(acc_fig3)
fig3 <- plot_ly(acc_fig3, x = ~MES, y = ~acc_fig3$`2018`, name = '2018',
type = 'scatter', mode = 'lines',
line = list(color = 'rgb(205, 12, 24)', width = 4))
fig3 <- fig3 %>% add_trace(y = ~acc_fig3$`2017`, name = '2017', line = list(color = 'rgb(22, 96, 167)', width = 4))
fig3 <- fig3 %>% add_trace(y = ~acc_fig3$`2016`, name = '2016', line = list(color = 'rgb(205, 12, 24)', width = 4, dash = 'dash'))
fig3 <- fig3 %>% add_trace(y = ~acc_fig3$`2015`, name = '2015', line = list(color = 'rgb(22, 96, 167)', width = 4, dash = 'dot'))
fig3 <- fig3 %>% add_trace(y = ~acc_fig3$`2014`, name = '2014', line = list(color = 'rgb(205, 12, 24)', width = 4, dash = 'dot'))
fig3 <- fig3 %>% layout(title = "Accidentes por Mes",
xaxis = list(title = "Meses"),
yaxis = list (title = "Accidentes"))
fig3
A partir de la serie de tiempo graficada se pueden comentar las siguientes ideas: 1. Los cuatro años concuerdan en que el mes con menos accidentes es Enero. Es un mes “de reposo” despues de las festividades de Fin de Año. 2. Los accidentes en el mes de marzo fueron disminuyendo año tras año. la direfencia es considerable al pasar en 2014 de 3855 accidentes a 3335 accidentes el 2018. 3.Se presenta el mismo pico en el mes Mayo en 3 de los 4 años; siendo por año, uno de los meses con mayor accidentalidad. 4. 3 de los 4 años, concuerdan en que Agosto es el mes con mayor accidentalidad del año.
Estos comportamientos pueden darse dado fechas especiales de cada mes con alta frecuencia de accidentes; por ejemplo en marzo el día de la mujer; en mayo el día de la madre, en agosto la feria de flores.
Con base en lo anterior, se hará una ampliación sobre la siguiente medida de tiempo, los días por mes. Con base en lo anterior, se hará una ampliación sobre la siguiente medida de tiempo, los días por mes. De esta mandera el siguiente mapa de calor mostrará cuáles son los días por mes que mayor accidentalidad presentan.
acc_fig4 <- as.data.frame(as.Date(accidentes[, "FECHA"], format = "%Y-%m-%d"))
colnames(acc_fig4) <- "FECHA"
acc_fig4$MES <- strftime(acc_fig4$FECHA, format = "%B")
acc_fig4$MES <- factor(acc_fig4$MES, levels = c("enero", "febrero", "marzo", "abril", "mayo", "junio", "julio", "agosto", "septiembre", "octubre", "noviembre", "diciembre"))
acc_fig4$DIA <- strftime(acc_fig4$FECHA, format = "%A")
acc_fig4$DIA <- factor(acc_fig4$DIA, levels = c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
acc_fig4m <- as.matrix(table(acc_fig4$MES, acc_fig4$DIA))
acc_fig4m <- apply(acc_fig4m, 2, function(x){x/mean(x)})
fig5 <- plot_ly(x=colnames(acc_fig4m), y=rownames(acc_fig4m),
showscale = F,
z = acc_fig4m,
colors = colorRamp(c("lightgoldenrod1", "indianred2")),
type = "heatmap")
fig5
La gráfica de calor muestra que los días viernes y jueves en los meses agosto y septiembre son los días con mayor accidentalidad en los 4 años considerados y que los días con menos accidentalidad son los lunes, específicamente en los meses enero, junio y noviembre; sin embargo el día que por mes tiene más accidentalidad es el sábado.
Habiendo visto de manera gráfica como se comportan las frecuencias de los accidentes, se procede a depurar y limpiar la base de datos con el objetivo de encontrar nuevas variables entre las relaciones de las existentes; además se continúa con la sección de agrupamiento.
## Agrupamiento
table(accidentes$GRAVEDAD)
##
## HERIDO MUERTO SOLO DAÑOS
## 114201 1264 93973
table(accidentes$CLASE)
##
## Atropello Caida ocupante Choque
## 6 20679 18578 141447
## Choque y Atropello Incendio Otro Volcamiento
## 1 24 22066 6637
Dado que la clase del accidente puede proveer estadísticas importantes para realizar el agrupamiento de barrios se procede a unificar para tener categorías con una cantidad importante de registros. Por esta razón, las categorías vacías e Incendio serán agregadas a la categoría Otro, por otro lado, La categoría Choque y Atropello se agregará a los datos de la categoría Choque.
accidentes <- accidentes %>%
rename(BARRIO_1 = BARRIO)
accidentes$CLASE <- gsub("Incendio","Otro",accidentes$CLASE)
accidentes$CLASE <- gsub("Choque y Atropello","Choque",accidentes$CLASE)
accidentes$CLASE[accidentes$CLASE == ""] <- "Otro"
table(accidentes$CLASE)
##
## Atropello Caida ocupante Choque Otro Volcamiento
## 20679 18578 141448 22096 6637
Con esta información es necesario idear algunos indicadores que más allá de sus valores aporten información y sean útiles para la intervención de la accidentalidad. En un principio se pueden establecer indicadores para cada una de las categorías de la clase de accidente así como la gravedad de los mismos y la población de cada barrio.
Hasta el momento los indicadores definidos para el agrupamiento son:
Ya que existen algunos datos faltantes para el tema de los barrios se procede a asignarle su respectivo barrio de acuerdo a sus coordenadas. Para se importa el archivo geoJson. que contiene los polígonos de los barrios de Medellín. Este procedimiento es igualmente válido para poder graficar el mapa adecuadamente con los resultados del clustering.
Sys.setlocale(locale = "Spanish")
## [1] "LC_COLLATE=Spanish_Spain.1252;LC_CTYPE=Spanish_Spain.1252;LC_MONETARY=Spanish_Spain.1252;LC_NUMERIC=C;LC_TIME=Spanish_Spain.1252"
barrios_medellin <- readOGR("G:/Mi unidad/UNIVERSIDAD/ESPECIALIZACION/Analítica Predictiva/TRABAJO/analitica_predictiva/Agrupamiento/Limite_Barrio_Vereda_Catastral.geojson",
stringsAsFactors = FALSE,use_iconv = TRUE,encoding = "UTF-8")
## OGR data source with driver: GeoJSON
## Source: "G:\Mi unidad\UNIVERSIDAD\ESPECIALIZACION\Analítica Predictiva\TRABAJO\analitica_predictiva\Agrupamiento\Limite_Barrio_Vereda_Catastral.geojson", layer: "Limite_Barrio_Vereda_Catastral"
## with 354 features
## It has 8 fields
Para realizar este proceso es necesario que los puntos provenientes de la accidentalidad se encuentren en la misma proyección de los polígonos.
accidentes_spatial <- accidentes
coordinates(accidentes_spatial) <- ~LONGITUD+LATITUD
proj4string(accidentes_spatial) <- proj4string(barrios_medellin)
## Warning in proj4string(barrios_medellin): CRS object has comment, which is lost
## in output
Ahora se llama la función over de la librería sp para determinar los poligonos en los que se encuentra cada punto, este proceso ayudará a declarar los nombres de los barrios de una manera estándar para generar el mapa.
## Código para determinar en cuáles barrios se encuentran los puntos
contains_barrio <- over(accidentes_spatial, barrios_medellin)
De esta manera, ahora se añade la información del archivo geojson al data frame para poder identificar los puntos que se encuentran sin barrio.
accidentes <- cbind(accidentes,contains_barrio)
Se someterá el set de datos a una búsqueda de datos faltantes:
table(is.na(accidentes$NOMBRE_BARRIO))
##
## FALSE TRUE
## 208923 515
Siguen existiendo 515 registros sin un barrio asignado, los cuales se apartan del set de datos base.
accidentes_na <- accidentes[is.na(accidentes$NOMBRE_BARRIO),]
Después de analizar los datos faltantes, se puede observar que estos registros se encuentran por fuera de la región considera Medellín de acuerdo a las delimitaciones encontradas en el archivo geojson, por esta razón se eliminarán del conjunto de datos.
accidentes <- accidentes[!is.na(accidentes$NOMBRE_BARRIO),]
table(is.na(accidentes$NOMBRE_BARRIO))
##
## FALSE
## 208923
Ahora con estos datos podemos proceder a realizar mapas y hacer el análisis por clustering. Para comenzar se ilustrará un ejemplo simple con el número total de accidentes.
## Obtener indicador de cantidad de accidentes por barrio en total
barrios_ejemplo <- accidentes %>%
select(NOMBRE_BARRIO) %>%
group_by(NOMBRE_BARRIO) %>%
summarize(conteo = n())
## `summarise()` ungrouping output (override with `.groups` argument)
head(barrios_ejemplo)
## # A tibble: 6 x 2
## NOMBRE_BARRIO conteo
## <chr> <int>
## 1 Aguas Frias 63
## 2 Aldea Pablo VI 36
## 3 Alejandría 428
## 4 Alejandro Echavarría 773
## 5 Alfonso López 1239
## 6 Altamira 720
## Copia del spatial data frame
barrios_med_ejemplo <- barrios_medellin
## Sacar index para organizar luego los registros
barrios_med_ejemplo@data$index <- as.integer(row.names(barrios_med_ejemplo@data))
## Combinar el indicador con el data frame spatial para graficar en mapa
barrios_med_ejemplo@data <- merge(barrios_med_ejemplo@data,barrios_ejemplo,by = "NOMBRE_BARRIO",all.x=TRUE)
barrios_med_ejemplo@data <- barrios_med_ejemplo@data %>% arrange(index)
barrios_med_ejemplo@data$conteo[is.na(barrios_med_ejemplo@data$conteo)] <- 0
Mapa
labels <- sprintf(
"<strong>%s</strong><br/>%g accidentes en total",
barrios_med_ejemplo@data$NOMBRE_BARRIO, barrios_med_ejemplo@data$conteo
) %>% lapply(htmltools::HTML)
pal <- colorNumeric("YlOrRd", NULL)
leaflet(barrios_med_ejemplo)%>%
addTiles() %>%
setView(lat=6.247612, lng=-75.582932, zoom=11.5) %>%
addPolygons(color="grey",weight = 1, fillOpacity = 0.75,opacity = 1, smoothFactor = 0.7, fillColor = ~pal(conteo),
highlight = highlightOptions(
weight = 3,
color = "#666",
fillOpacity = 0.8,
bringToFront = TRUE),
label = labels,
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "15px",
direction = "auto")) %>%
addLegend(pal = pal, values = ~conteo, opacity = 0.7, title = "Cant. accidentes",
position = "topright")